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The structural properties of the ZnxMgi_xSySei_y sohd solutions are determined by a combination 
of the computational alchemy and the cluster expansion methods with Monte Carlo simulations. We 
determine the phase diagram of the alloy and show that the homogeneous phase is characterized by 
a large amount of short-range order occurring among first-nearest neighbors. Electronic-structure 
calculations performed using the special quasi-random structures approach indicate that the energy 
gap of the alloy is rather sensitive to this short-range order. 



Wide gap semiconductors have recently attracted an 
enormous technological interest 0| both because of their 
potential use in devices capable of operating at high 
power level and high temperature, and because of the 
need for optical materials active in the blue-green spec- 
tral range. ZnSe-based technology will be used for oper- 
ation in this spectral range provided that current device 
lifetime problems are overcome. One major goal of ma- 
terials engineering for opto-electronic applications is the 
ability to tune independently the band gap, Eg — in order 
to obtain the desired optical properties — and the lattice 
parameter, oq, of the material — in order to be able to 
grow it on a given substrate. Unfortunately, in most 
III-V and II- VI alloys, the additional degree of freedom 
provided by alloying both the cationic and the anionic 
sub-lattices cannot be effectively exploited because Eg 
and flp are roughly inversely proportional to one another 
for any values of the cationic and anionic compositions, 
(x,2/). From this point of view, ZnxMgi_xSySei_y alloys 
play a special role, in that the lattice parameter and op- 
tical gap can be varied fairly independently as functions 
oiix,y)_§. 

In spite of this big interest, many technical difficul- 
ties still hinder a precise experimental characterization 
of these materials, so that their equilibrium structural 
and optical properties are basically unknown. In this 
Letter we report on the first application of state-of-the- 
art electronic-structure techniques to the determination 
of the structural and optical properties of a quaternary 
(double binary) semiconductor alloy at thermodynamic 
equilibrium, and present results in the specific case of 
ZnxMgi_xSySei_y. Even though MBE-grown materials 
are fabricated in highly non-equilibrium conditions, un- 
derstanding their equilibrium properties is preliminary 
to any further investigations and provides the physical 
limits to the tunability of the alloy properties. 

The first goal of this Letter is to determine the ther- 
modynamic stability of the (Zn,Mg)(S,Se) solid solution 
with respect to segregation into its constituents and/or 
to the formation of ordered structures. Secondly, we will 
analyze how the fundamental gap depends on composi- 
tions and the role that short-range order plays in the 



electronic properties of this material. 

The thermodynamic properties of ZnxMgi_xSySei_y 
are studied by mapping the alloy onto a (double) lattice- 
gas model [HQ that is solved by standard Monte Carlo 
techniques. To this end, an Ising-like variable, {ctrs}, 
is first attached to the s-th atom in the elementary cell 
located at R, and it is assumed to take the values ±1 
according to the type of atom occupying that lattice site. 
The energy of the alloy is then expanded in terms of the 
(t's including terms up to three-spin interactions: 
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The validity of the truncation is then established a poste- 
riori. The two- and three-body interaction constants, J 
and L, can in principle be obtained from both the cluster 
expansion (CE) [^ or computational alchemy (CA) jj^] ap- 
proaches. The explicit inclusion of local distortions from 
the ideal lattice positions — due to the different atomic 
sizes — renormalizes the two-body interaction, making it 
long-ranged [g[ and hence hardly obtainable by the CE 
method. In principle, the three-body interactions, L, 
could also be computed within CA by making use of the 
so called 2n+l theorem |j]. Assuming that the L's are 
short-ranged, however, they are more conveniently ob- 
tained by the CE method. To this end, we first calculate 
the J's up to 9-th nearest neighbors from CA by using 
second-order density- functional perturbation theory [pj. 
We then assume that the only non- vanishing L's are those 
which connect atoms that are first- or, at most, second- 
nearest neighbors. Because of symmetry, there are only 
six such 3-body interaction constants. In the zincblende 
structure, there are two inequivalent triplets of second- 
nearest neighbors which differ according to whether or 
not they have a common first-nearest neighbor. We fol- 
low the common practice of neglecting such difference 
[H, and we keep 4 inequivalent 3-body interactions. In 
order to determine the latter, we have calculated the total 



energies of a set of 28 pseudo-binary ordered structures 
by using self-consistent density-functional theory (DFT) 
within the local-density approximation (LDA) [pi, and 
by allowing the lattice to relax until the forces acting on 
individual atoms vanish. We have then fitted the cubic 
term of Eq. (1) to the differences between the energies 
so calculated and the predictions of Eq. (1), when trun- 
cated to second order. The quality of the fit so obtained 
is finally checked against 39 additional quaternary struc- 
tures. The resulting mean-square error is ^ 0.8 meV per 
atom pair, while the maximum error is 2.1 meV. Two- 
and three-body constants have been calculated for 5 and 
3 volumes in the range [T^znSi ^MgSc] respectively, and in- 
terpolated in-between. 

Finite-temperature properties have been determined 
by lattice-gas Monte Carlo simulations, using supercells 
of 1024 atoms, at constant temperature, pressure, and 
chemical potentials. To this end, we attempt two dif- 
ferent types of Monte Carlo moves: the reversal (flip) 
of each spin and a variation of the volume of the super- 
cell. The trial move is then accepted with the Metropolis 
probability. Measures were taken over « 10'^ correlation 
times after thermal equilibrium was reached. 

The free energy surface of the system is extracted from 
the determination of the compositions, x and y, as func- 
tions of the chemical potentials, fix and /i^. The regions 
of spinodal (local) stability correspond to those concen- 
trations for which the Hessian of the free energy with re- 
spect to the concentrations is positive definite. The ther- 
modynamically stable concentrations (binodal regions), 
and thus the miscibility gaps, are determined with a gen- 
eralization of the common-tangent Maxwell construction 
usually adopted for binary mixtures: the globally stable 
points are those which are locally stable and whose tan- 
gent planes do not intersect the free energy surface at 
any other points in the square of concentrations. 

In Fig. |l| we display a cut of the phase diagram of the 
quaternary alloy at four different temperatures. At typi- 
cal MBE growth temperatures (« 600 K) and above the 
alloy is predicted to be completely miscible. For T = 
550 K (Fig. pi), the disordered alloy is stable at all com- 
positions but for a small island close to the ZnxMgi_xSe 
pseudo-binary alloy centered around x w 0.6, whose com- 
puted critical temperature is of about 610 K. Decreasing 
the temperature (Fig. |^b) the size of the island increases 
and a forbidden region appears inside the square of com- 
positions, close to the mid point {x,y) = (ji 2)' Cooling 
down to T = 500 K (Fig. Qc), the Mg-rich and the Zu- 
rich regions become separated by a "corridor" in which 
the disordered alloy is not locally stable. The concen- 
trations of the phases in which the alloy segregates in 
the miscibility gap are given by the contact points of the 
tangent plane. Segregation may thus result in the sep- 
aration into two or three phases, according to the num- 
ber of contact points. In the present case, the critical 
temperatures of the pseudo-binary alloys with cationic 
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FIG. 1. Phase diagram of the alloy at four different tem- 
peratures. The dark gray (binodal) regions are thermody- 
namically stable. The light gray (spinodal) regions are only 
locally stable (metastable). The white regions are completely 
unstable. The segments inside the miscibility gap connect the 
stable phases identified by the common-tangent-plane con- 
struction (see text). In panels a — c only separation into two 
phases occur. In panel d two cases where the products of 
separation are three distinct phases are indicated by thicker 
lines. 

disorder, ZnxMgi_xSe and ZnxMgi_xS, are much larger 
{Tc = 613 K and 511 K respectively) than those of the 
pseudo-binary alloys with anionic disorder, ZnSySei_y 
(Tc = 254 K) and MgSySei_y (T^ = 243 K); therefore, 
above room temperature the alloy segregates into two 
phases: a Zn-rich phase and a Mg-rich one. This be- 
havior is due to the larger chemical difference between 
cations than between anions. At temperatures of the or- 
der of the anionic-alloy critical temperatures and below 
(T < 250 K), the tendency to segregation also involves 
the anionic sub-lattice, and separation into three phases 
may occur (Fig. |jd). 

Our results show that in (Zn,Mg)(S,Se) the onset of 
segregation occurs at a critical temperature that is in 
the range of typical MBE growth temperatures, and thus 
much lower than in other II- VI quaternary alloys f^. 
This result is in agreement with an analysis of experi- 
mental data [ [lO[ , based on delta-lattice-parameter mod- 
els [0, that locate T^ between 525 K and 625 K. Ac- 
cording to the only experimental report p2[ we are aware 
of, a forbidden region of compositions has been observed, 
at room temperature, inside the predicted miscibility 
gap. 

The structural properties of the alloy can be extracted 
from our Monte Carlo simulations. The equilibrium lat- 
tice constant has a strong linear dependence upon com- 
positions, thus following Vegard's law. Slight deviations 
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LDA 

Corrected LDA 



11.326 
11.130 



-0.789 
-0.443 



-0.539 
-0.523 



0.058 
0.018 



0.027 
0.027 



0.016 
0.016 



TABLE I. Quadratic fit coefficients for LDA (upper row) 
and "corrected" LDA (lower row) lattice constant of the alloy. 
The maximum error of the fit is of 0.05%, and the mean square 
error of 0.01%. 
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FIG. 2. Left panel: bond-length distributions along the 
ZnS — MgSe diagonal in the square of the compositions. Right 
panel: dependence of maxima of the peak upon the concentra- 
tions. Diamonds correspond to pure-material bond-lengths. 
Empty dots refer to the alloy lattice constant and the dashed 
line is a linear fit. 

from this law are extracted by a fit of the lattice param- 
eter of the alloy with a quadratic polynomial in x and 
y (see Table ffl). A direct comparison with experimental 
data may be done by correcting for the LDA error on the 
equilibrium lattice constants of the pure constituents of 
the alloy (±1 ~ 2%), which is accounted for by a bilin- 
ear interpolation and added to the quadratic polynomial 
obtained above. The corrected coefficients of polynomial 
fit are also reported in Table ffl. The agreement with 
experimental data is good in the Zn-rich region, es- 
pecially in proximity of pure ZnSe, while it is apparently 
poor in the Mg-rich part of the square of the composi- 
tions, where however the instability of the system makes 
it very difficult to measure the concentrations with the 
same precision, and the experimental data are in fact 
very scarce. 

As it is the case in other tetrahedrally bonded solid so- 
lutions, the Vegard's law does not result from a linear de- 
pendence of individual bond lengths upon compositions, 
but rather from a subtle topological compensation of 
bond lengths which stay in fact rather close to the values 
they would have in pure compounds Eol. This behavior 
is displayed in Fig. 0, where we show the bond lengths as 
obtained at T=800 K for the (Zn,Mg)(S,Se) alloy along 
the MgSe-ZnS diagonal of the square of concentrations 
(along the other diagonal the lattice parameter is almost 
constant and matched to that of GaAs). At any concen- 
trations, while the lattice parameter of the alloy varies 
almost linearly by more than 12 %, the largest deviations 
between the bond lengths and their pure-compound val- 
ues are smaller than 3%. 
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FIG. 3. First- and second-nearest-neighbors two-body cor- 
relations of the Zn i Mg i S i Se i alloy as function of tempera- 
ture. 

The typical substrate on which such alloys are grown 
is GaAs, whose lattice constant is ao — 10.68 a.u., close 
to the equilibrium lattice parameters of ZnSe and MgS. 
Therefore, the technological relevant concentrations of 
the quaternary alloy are those located along the ZnSe- 
MgS diagonal of the square of the compositions. For 
these solid solutions, the Mg-Se and Zn-S bonds are con- 
sequently subject to a larger elastic strain with respect 
to the other bonds; the preferential formation of MgS 
and ZnSe clusters is thus expected to be energetically fa- 
vored. In order to clarify this issue we have calculated 
the two-body correlations of the ZniMgiS i Sei alloy at 
temperatures above the miscibility gap. The correlation 
function is defined as: Css'(R) = (ctsoCs'r) — {(Ts){<J8')- 
Therefore, Css'(R) = in the perfectly random alloy, 
Css'(R.) > for ZnSe and MgS clusterizations, and 
Css' (R.) < for ZnS and MgSe clusterizations. In corre- 
spondence of the first-nearest-neighbor shell, Css' (R) has 
a very pronounced positive peak, thus indicating a strong 
tendency to form ZnSe and MgS clusters. This tendency 
is confirmed by the positive value of the cation-cation or 
anion-anion correlations in the second-nearest-neighbor 
shell, which are however weaker. Correlations practically 
die beyond the second shell of neighbors, and are there- 
fore characteristic of a very pronounced short-range or- 
der (SRO). In Fig. ^we display the value of the first- and 
second-shell correlation peaks for a range of temperatures 
above Tc- The nearest-neighbor correlation peak is still 
significantly large at very high temperature, thus indi- 
cating that short-range-order is present even close to the 
melting temperatures (above T k. 1700 K): the system 
can never be described as a perfectly random alloy. 

The second goal of this work is to study the influence of 
the structural properties of the alloy at the atomistic level 
on the electronic and optical properties. The band struc- 
ture of a disordered material is strongly affected by its lo- 
cal environment that is poorly approximated by effective- 
medium approaches |l^,|lj]. A detailed account of the 
microscopic structure of the alloy is even more impor- 
tant in the present case where SRO effects are expected 
to play a significant role. A "direct" DFT-LDA calcula- 
tion should be in principle performed by using very large 
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FIG. 4. Energy gap of the (Zn, Mg)(S, Se) quaternary 
alloy lattice-matched to GaAs as function of the x composi- 
tion (bottom scale) and y composition (top scale). The dot- 
dashed line refers to virtual crystal calculations, the dashed 
one to perfectly-random-alloy results, and the solid line to the 
short-range-ordered alloy band gaps. 

supercells in order to cope with compositional disorder. 
Here we adopt the special quasi-random structures (SQS) 
approach |^ , suitably generalized to account for SRO ef- 
fects [E6J, double-sublattice disorder, and arbitrary com- 
positions. We have considered three different concentra- 
tions along the line of lattice matching to GaAs, and 
chosen among the discrete set of compositions compati- 
ble with 64-atom SQS's and close to the ZnSe-rich region: 
{x,y) = (^, 5); (|, |); (HjO). For these concentrations, 
SQS's have been obtained by a simulated-annealing pro- 
cedure aimed at modeling the pair correlation functions 
at short range — as obtained from Monte Carlo simula- 
tions performed at T = 550 K and resulting in very 
faithful pair correlations up to fourth-nearest-neighbors. 
We have explicitly verified that the calculated electronic 
properties of the alloy are rather insensitive to correla- 
tions beyond this order of neighbors. 

It is interesting to compare the results obtained using 
different levels of approximations to substitional disor- 
der. In Fig. we display the energy gap as a function of 
the cationic composition along the GaAs-matching line 
of the compositions plane, as obtained from calculations 
performed i) on the appropriate virtual crystal, ii) on 
supercells describing the perfectly random alloy, and Hi) 
on SQS's that reproduce the SRO correlations as ex- 
plained above. We see that virtual-crystal is a bad ap- 
proximation of the real system, while the effects of SRO 
show up in a slight but non-negligible opening of the 
fundamental band gap. Our results, once bilinearly cor- 
rected for the DFT error in the pure materials, agree very 
well (within 40 meV), with freshly appeared experimen- 
tal estimates |l7|. An analysis of the momentum- and 
position-projected density of states shows that the fun- 
damental gap is direct and occurs at the F point of the 
Brillouin zone for any concentrations. The top-valence 
band states are mainly Se states, while the bottom state 
of the conduction band is rather delocalized on the differ- 
ent atomic species and it is found to be the most sensitive 



to the occurrence of short-range-order. 

We wish to thank A. Franciosi and S. Rubini for 
prompting our interest in this problem and for many 
fruitful discussions. A more complete account of this 
work can be found in the PhD thesis of one of us (AMS) 
which is available on the www |1n] . This work has been 
done in part within the Iniziativa Trasversale Calcolo 
Parallelo of the INFM. 

" Present address: Center for Molecular Modeling, Dept. 

of Chemistry, University of Pennsylvania, 231 S 34th St, 

Philadelphia, PA 19104-6202 
* Electronic addresses: saitta@cmm.chem.upenn.edu, de- 

gironc@sissa.it, and baroni@sissa.it. 
[1] II-VI Blue/Green Light Emitters: Device Physics and 

Epitaxial Crowth, ed. by R. L. Gunshor and A. V. Nur- 

mikko. Semiconductor and Semimetals 44, Academic 

Press, (1997). 
[2] H. Okuyama, K. Nakano, T. Miyajima, and K. Akimoto, 

Jpn. J. Appl. Phys. 30, L1620 (1991); J. Cryst. Growth 

117, 139 (1992). 
[3] D. de Fontaine, in Solid State Physics, edited by H. 

Ehrenreich, F. Seitz, and D. TurnbuU (Academic, New 

York, 1979), 34, 73. 
[4] J. W. D. Connolly and A. R. Wilhams, Phys. Rev. B 27, 

5169 (1983). 
[5] L. G. Ferreira, S.-H. Wei, and A. Zunger, Phys. Rev. B 

40, 3197 (1989); S.-H. Wei, L. G. Ferreira, and A. Zunger, 

Phys. Rev. B 41, 8240 (1990). 
[6] S. de Gironcoli, P. Giannozzi, and S. Baroni, Phys. Rev. 

Lett. 66, 2116 (1991); N. Marzari, S. de Gironcoh, and 

S. Baroni, Phys. Rev. Lett. 72, 4001 (1994). 
[7] X. Gonze and J. P. Vigneron, Phys. Rev. B 39, 13120 

(1989). 
[8] All our ab-initio calculations have been performed in the 

plane-wave pseudopotential scheme, adopting the non- 
linear core-correction description for Zn and Mg atoms, 

and choosing a kinetic-energy cutoff of 15 Ry and special 

k-points sampling of the Brillouin zone. 
[9] D. W. Kisker, J. Cryst. Growth 98, 127 (1989). 
[10] T. Maruyama, et al, J. Cryst. Growth 159, 41 (1996). 
[11] G. B. Stringfellow, J. Cryst. Growth 27, 21 (1974); ibid. 

65, 454 (1983). 
[12] B. J. Wu, et al, Appl. Phys. Lett. 66, 3462 (1995). 
[13] L. Nordheim, Ann. Phys. (Leipzig) 9, 607 (1931). 
[14] P. Soven, Phys. Rev. 156, 809 (1967); D. W. Taylor, ibid. 

156, 1017 (1967); U. Onodera and Y. Toyozawa, J. Phys. 

Soc. Jpn. 24, 341 (1968); B. Velicky, S. Kirkpatrick, and 

H. Ehrenreich, Phys. Rev. 175, 747 (1968). 
[15] A. Zunger, S.-H. Wei, L. G. Ferreira, and J. E. Bernard, 

Phys. Rev. Lett. 65, 353 (1990); S.-H. Wei, L. G. Ferreira, 

J. E. Bernard, and A. Zunger, Phys. Rev. B 42, 9622 

(1990); K. C. Hass, L. C. Davis, and A. Zunger, Phys. 

Rev. B 42, 3757 (1990). 
[16] K. A. Mader and A. Zunger, Phys. Rev. B 51, 10462 

(1995). 
[17] U. Lunz, et al, Semicond. Sci. Technol. 12, 970 (1997). 
[18] A. M. Saitta, SISSA Ph. D. thesis, available at the URL: 



tittp: //www. sissa. it/cm/thesis/1997/ 



